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<^ i Abstract 



Recent experiments have shown that the superconducting energy gap in some cuprates is spatially 

o 

O ' inhomogeneous. Motivated by these experiments, and using exact diagonalization of a model d- 

£4 

Q_i' wave Hamiltonian, combined with Monte Carlo simulations of a Ginzburg-Landau free energy 

2. 

functional, we have calculated the single-particle density of states LDOS(cj,r) of a model high-T c 

+-> 

superconductor as a function of temperature. Our calculations include both quenched disorder in 
the pairing potential and thermal fluctuations in both phase and amplitude of the superconducting 
gap. Most of our calculations assume two types of superconducting regions: a, with a small gap 

o : 

and large superfluid density, and f3, with the opposite. If the p regions are randomly embedded in 
an a host, the LDOS on the a sites still has a sharp coherence peak at T = 0, but the (3 component 
does not, in agreement with experiment. An ordered arrangement of (3 regions leads to oscillations 



> 

in 
in 
o . 

in the LDOS as a function of energy. The model leads to a superconducting transition temperature 



\Q ' T c well below the pseudogap temperature T c q, and has a spatially varying gap at very low T, both 

o ' 

consistent with experiments in underdoped Bi2212. Our calculated LDOS(oj,r) shows coherence 



a 

but not amplitude fluctuations in a homogeneous superconductor. Well above T c , the gap in the 

O 
O 

> 



peaks for T < T c , which disappear for T > T c , in agreement with previous work considering phase 
but not amplitude 
LDOS disappears. 
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I. INTRODUCTION 



According to low temperature scanning tunneling microscopy (STM) experiments, 
;he local density of states (LDOS) of some cuprate materials have spatial variations 



flflflflflflflg. 



Among the cuprates, E^S^CaC^Og+x (Bi2212) is one of the 
most extensively studied in STM experiments. The LDOS spectrum shows that some re- 
gions of that material, which we will call a-regions, have a small energy gap with large 
and narrow coherence peaks (reminiscent of the spectra observed in bulk superconducting 
materials), while other regions, which we will call /3-regions, have a larger gap, but smaller 
and broadened peaks (which are reminiscent of the spectra seen in bulk pseudogap phase 
of some materials . These inhomogeneities occur on length scales of order 30A. Because at 
low doping concentrations a regions with "good" superconductivity are immersed in more 
metallic or semiconducting (3 regions, some workers have made an analogy between these 
materials and granular superconductors 0,0]: superconducting domains spatially separated 
from one another by non-superconducting regions, but connected through proximity effect 
or Josephson tunneling. 

At present there is no general agreement regarding the origin of the inhomogeneities in 
the cuprate superconductors — whether they are in charge density, spin density, LDOS, or 



other properties |l0|. One hypothesis is that these in 



self organization due to competing orders 



romogeneities originate in a process of 



111, 12, 13, 14, Q 16, Ell- 



in another approach, 



the spatially varying properties of the cuprates are attributed to crystal defects or impurities. 
In particular, it has been suggested that the inhomogeneities in the LDOS originate in the 



random spatial distribution of dopant atoms near the copper oxide (Cu0 2 ) planes 
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18, 



2l|. 



Several workers have studied the LDOS of inhomogeneous superconductors at low T. 
For example, Ghosal et al{22] have calculated the LDOS of a strongly disordered s-wave 
superconducting layer in two dimensions, solving the Bogoliubov-de Gennes equations self- 
consistently. They have also done similar work on a model of c?-wave superconductivity 
Fang et a/.Q, using a Green's function approach, computed the zero temperature LDOS 
of a model lattice Hamiltonian in which one small region of the lattice has an different 
(either suppressed or enhanced) pai ring strength than the rest; they find good agreement 
with experiments. Cheng and Su |24| have also explored how the LDOS is affected by 
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a single spatial inhomogeneity in the pairing strength of a BCS Hamiltonian; they find 
that an inhomogeneity with an LDOS most closely resembling the experimental results is 
produced by an inhomogeneity with a cone-shaped distribution of the pairing strength; this 
work thus suggests that it is the small-length-scale variation of the pairing strength that 
causes incoherence in the LDOS. Mayr et al. |25( have studied a phenomenological model 
with quenched disorder and observed a pseudogap in the LDOS caused by a mixture of 
antiferromagnetism and superconductivity, while Jamei et al. j3] have investigated the low 
order moments of the LDOS and their relation to the local form of the Hamiltonian. 

In this paper we propose a phenomenological approach to study the effect of inhomo- 
geneities on the LDOS in a model for cuprate superconductors. The model is a mean-field 
BCS Hamiltonian with <i-wave symmetry, in which the pairing-field is inhomogeneous and 
also undergoes thermal fluctuations in both phase and amplitude at finite temperatures T. 
It has been argued jl(3] that the superconducting state of optimally doped to overdoped 
cuprates is well described by the BCS theory which includes a d-wave gap and scatter- 
ing from defects outside the Cu02 plane. Instead of including such defects explicitly in our 
BCS Hamiltonian, we implicitly include their possible effects through inhomogeneities of the 
pairing-field amplitude. Furthermore, instead of self-consistently solving the Bogoliubov-de 
Gennes equations resulting from this model, we obtain the magnitude and phase of the com- 
plex pairing-field from Monte Carlo simulations based on a Ginzburg Landau free energy 
functional. Thus the procedure is as follows. First, we set the parameters of the Ginzburg 
Landau free energy functional from experiments. Next, using Monte Carlo simulations of 
this free energy, we obtain the pairing-field amplitudes which we then include in the BCS 
Hamiltonian. Finally, we diagonalize the latter in order to obtain the LDOS. 

Now in optimally or nearly optimally doped Bi2212, the layers consist of randomly dis- 
tributed /3-regions immersed in a majority background of a-regions (^J. We therefore choose 
Ginzburg-Landau parameters so as to reproduce this morphology at T = 0, then carry out 
simulations at both zero and finite T to obtain the LDOS in the different spatial regions. 

At T = we compare these simulation results to those obtained using ordered instead of 
random arrangements of inhomogeneities. We find that the LDOS of the random systems 
much more closely resemble experiment. Specifically, regions with a small gap have sharp 
coherence peaks, while large-gap regions show lower and broader peaks. By contrast, systems 
with ordered inhomogeneities have LDOS spectra with sharp coherence peaks which oscillate 
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as a function of energy. In the ordered systems, the coherence peaks in the small-gap regions 
strongly resemble those observed in a homogeneous small-gap system. But the spectral peaks 
in the large-gap regions dramatically differ from those in the corresponding homogeneous 
and disordered cases. 

Because the spectra of disordered systems more closely resemble experiments, we have 
also studied the evolution of the LDOS in these systems with increasing T. We consider 
both T < T c and T > T c , where T c is the phase-ordering transition temperature (equivalent 
to the Kosterlitz-Thouless transition temperature for this two-dimensional system). In both 
the a and j3 regions, we find that the spectral gap starts to fill in as T increases, and the 
spectral peaks broaden and are reduced in height. However, even above T c , the LDOS is 
still suppressed at low energies, in comparison to the normal state. This result agrees with 
a previous study j^i |3| which considered thermal fluctuations of the phase but not of the 
magnitude of the complex pairing-field, and included no quenched disorder. 

We have also studied the T-dependence of the magnitude of the pairing field, its thermal 
fluctuations, and the effective superfluid density of our disordered system. We find that the 
phase ordering temperature is greatly reduced from the spatial average of the mean-field 
transition temperatures appearing in the Ginzburg Landau free energy functional. This 
reduction is due to both thermal fluctuations and quenched disorder in our model. 

Although our work involves a non-self-consistent solution of a (i-wave BCS Hamiltonian, 
it differs from previous studies of this kind[7, 
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23, 241 because it includes thermal fluctu- 



ations as well as quenched disorder in the pairing-field amplitude. For our model, quenched 
disorder is crucial in obtaining LDOS spectra which depend smoothly on energy and are 
also consistent with the observed low and broad peaks in the (3 regions. 

The rest of this article is organized as follows: In Section II, we present the BCS model 
Hamiltonian. In Section III, we derive the discrete form of the Ginzburg Landau free energy 
functional used in our calculations. We also discuss simple estimates of the phase ordering 
temperature, our choice of model parameters and our method of introducing inhomogeneities 
into our model. Section IV describe the computational methods used at both zero and finite 
temperature. These methods include a classical Monte Carlo approach to treat thermal 
fluctuations, exact diagonalization to obtain the LDOS, and the reduction of finite size 
effects on the LDOS by the inclusion of a magnetic field. Section V presents our numerical 
results at both T = and finite T. A concluding discussion and summary are given in 
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Section VI. 



II. MODEL 



A. Microscopic Hamiltonian 

We consider the following Hamiltonian: 

Hbcs = 2 tj i c *W + 2 $^( A u c il c iT + c - c - V X! i ^ W 

Here, denotes a sum over distinct pairs of nearest neighbors on a square lattice with 

N sites, Cj a creates an electron with spin a (f or j) at site // is the chemical potential, Ay 
denotes the strength of the pairing interaction between sites i and j, and £y is the hopping 
energy, which we write as 

t%j I' hop- (2) 

where t hov > 0. 

n 

Following a similar approach to that of Eckl et al. |27[ , we take Ay to be given by 

1 I Ail + | A,- 1 . , 

*V = 7 o e °» ( 3 ) 



where 



and 



9j + %)/2, if bond (i, j) is in x-direction, 
% <( (4) 
9j + ^')/2 + 7T, if bond is in y-direction, 



A J = |A,|e^, (5) 



is the value of the complex superconducting order parameter at site j. We will refer to 
the lattice over which the sums in $IJ are carried out as the atomic lattice (in order to 
distinguish it from the XY lattice, which will be described in the next section.) The first 
term in Eq. (^J) thus corresponds to the kinetic energy, the second term is a BCS type of 
pairing interaction with <i-wave symmetry, and the third term is the energy associated with 
the chemical potential. 
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Eq. may also be written 



H BCS = * - Nfi, 



(6) 



where 



and 



* = 



z = 1,N 



-»1 



.4 



i A* 

A -t* 



(7) 



(8) 



Here £ and A are N x N matrices with elements £y [£y = £y, as given by eq. (J2J if i and j 
are nearest- neighbors, = —\i if z = j, and £y = otherwise] and Ay [Ay = Ay, as given 
by eq. (j3J) if i and j are nearest-neighbors, and Ay = otherwise]. 
Let U be the unitary matrix that diagonalizes A, i. e., 



We can then rewrite (EI) as 



with $ defined by 



B = U^A U, B diagonal. 



H BCS = $ - iV/i, 



(9) 
(10) 
(11) 



If we make the following definitions: 



$ = 



7iT 
7l 



z = l,iV 



(12) 



and 



U 



(13) 



Ujin) -v*{n) 

Vj{ri) u*{ri) 

where z labels the row and j the column oi N x N matrices, then we can see that (jllj) is the 
typical Bogoliubov - de Gennes transformation J^l^: 

TV 



3=1 
N 

3=1 



(14) 



Thus $ is a 2A-dimensional column matrix and U is a 2N x 2A-dimensional square matrix. 

Denoting the diagonal elements of the matrix B by E n , we can use (JHJ), © and (fT3*|) to 
obtain 



t A* 




u n {ri) 


= E n 


u n {ri) 






A -t* 




v n {n) 




v n (n) 



Eq. (JT5|) is the eigenvalue problem which must be solved in order to compute the local 
density of states, as we describe next. 



B. Explicit expression for the local density of states 

We wish to compute the local density of states, denoted LDOS(u;,rj), as a function of 
the energy uj and lattice position = (xi,yi) at both zero and finite temperature T. Given 
the value of the of the superconducting order parameter Aj at each lattice site, the matrix 
A can be constructed and the LDOS(co>,r, {Aj}) can be computed through j^] 

LDOS(cu,r ?; ,{A i })= \\u n {ri)\ 2 5{u} - E n ) + \v n {r i )\ 2 5{uo + E n )} (16) 

n,E n >0 

At T = all the phases 9{ are the same, since this choice minimizes the energy of the 
superconducting system. Thus, in this case, once we know {|Aj|} we can solve eq. (JT3J) for 
u n (ri), u n (r») and E n , and use this solution in (J16j) . At finite T, since Aj will thermally 
fluctuate, we need a procedure to obtain an average of ([16)1 over the relevant configurations 
of {A,}. We explain that procedure next. 



III. MODEL FOR THERMAL FLUCTUATIONS 



At finite T we compute LDOS(co>, r*) by performing an average of LDOS(co>, Tj, {Aj}) over 
different configurations {Aj}. Those configurations are obtained assuming that the thermal 
fluctuations of {Aj} are governed by a Ginzburg-Landau (GL) free energy functional F, 
which is treated as an effective classical Hamiltonian. 

The Ginzburg Landau free energy functional has been widely studied and applied to a 
variety of systems. It has been extensively used to study granular conventional supercon- 



ductors 
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34. 
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36 



37 



381 ] . Other studies have focused on the use of GL theory to 



describe the phase diagram of extreme type II superconductors 



391 ] . the influence of defects 
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411 ]. and the effect 



38, 



421- 



on the structure of the order parameter of <i-wave superconductors 
of thermal fluctuations on the heat capacity of high temperature superconductors 
Yet other researchers have derived the GL equations for vortex structures from microscopic 
theories [43]. There has also been interest studying the nature of the transition in certain 



parameter ranges for this type of model 



44, 



45]. 



In this section we discuss a procedure for obtaining a suitably discrete form of F, and 
determining its coefficients from experiments. [The final form of F is given by eq. ([29|h] We 
also discuss a way to estimate the phase ordering temperature using this model, the choice 
of the parameters that determine the GL coefficients, and finally a method of introducing 
inhomogeneities into the model. 



A. Discrete form of the Ginzburg-Landau free energy 

For a continuous superconductor in the absence of a vector potential, the Ginzburg- 
Landau free energy density has the form 

F . =«(£-i) i^+^r+J-^f. (it) 

Since \ip'\ 2 and F' have dimensions of inverse volume, and energy per unit volume, it follows 
that a and b have dimensions of energy, and (energy x volume), respectively. 

The squared penetration depth A 2 (T) and zero-temperature Ginzburg-Landau coherence 
length £ are related to the coefficients of F by j^] 

h 2 

a = — 
2 

and 



a = > ( 18 ) 



b = 8^1 (M) 2 (19) 

where p? B ~ 5.4 x 10 _5 eV-A 3 is the square of the Bohr magneton. 

Let us assume that the position-dependent superconducting energy gap A$ at is related 
to iph as in conventional BCS theory, through 

a,: 



r\2 



A,; 2 



(20) 



9.386; 

where we have also assumed that T c0 , b, and a are functions of position. The validity of (|2U|) 
can be verified by noting that in the absence of fluctuations F' is minimized by 



cOi 



Combining flU} with flUJ), we obtain, at T = 0, 

|A,(0)| 2 = 9.38(A; B T c0i ) 2 . 



(22) 



This result agrees well with experiment provided (i) T c q is interpreted as the temperature 
at which an energy gap opens according to ARPES experiments, and (ii) A(0) is taken as 
the low-temperature (T << T c ) magnitude of the gap observed in ARPES and tunneling 
experiments \4:6^ . 

In order to obtain a discrete version of the free energy functional, we integrate the free 
energy density (|17|) over volume to yield the free energy 



F'dV. 



(23) 



Assuming that ip' ~ constant within a volume ^d (where £o is the zero-temperature coher- 
ence length and d is thickness of the superconducting layer), we can discretize the layer into 
M cells of volume t^d. Using (J2~U|) . we can then write 



M 



- = El — 



A 2 (0) 



A,- 



M 



E 



1 



1 



+ 



t 2(9.38) A?(0) 
A, : 



A; 



kBT c oi 
A, 



Xi(0)k B T c0i \j(0)k B T c 



cOj 



where 



tfd 



(24) 



(25) 



32(9.38)vrm* 2 /i| 

l K\ ~ 2866 eV-A 2 if d = 10A. Except for d, K\ is independent of material-specific parame- 
ters. 

In ()24j) the sums are performed over what we will call the XY lattice, which is not 
necessarily the same as the atomic lattice used in ([TJ. In (J24|) . A,; = | A. i |e * is the value 
of the superconducting order parameter on the ith cell of the XY lattice. The third sum is 
carried out over distinct pairs of nearest-neighbors cells (ij). 

In order to see how the XFlattice and the atomic lattice are related, we now analyze some 
of the relevant length scales in our problem. Typically, the linear dimension of the XY lattice 
cell is taken to be the T = coherence length £o of the material in the superconducting 
layer. In a cuprate superconductor, e. g., Bi2212, £ ~ 15 A, while the lattice constant of 
the microscopic (atomic) Hamiltonian of eq. (1) - i. e., the distance between the Cu sites in 



the C11O2 plane - is ao ~ 5AA |46j. Thus, in this case, a single XY cell would contain about 
nine sites of the atomic lattice, on each of which the superconducting order parameter would 
have the same value A*. 

It is convenient to introduce a dimensionless superconducting gap 



^0 



and a dimensionless temperature 



t = 



k B T 



E 



(26) 



(27) 



where £0 is an arbitrary energy scale which will be specified below. We can then rewrite 

as 

M / . \ , W ! 

^(°)*c0i 



F A / t \ 1 2 1_ 

^ - ^ U* J A?(0)£k M + 2-, 2(9 .3 8) A 2 (0) - 



.4 IV>< 



Xi(0)t. 



cQi 



i>3 



A,(0)t 



cOj 



Aj(0)t c ojA : ,(0)t C Qj 



cos 



(28) 



In our calculations, we will employ periodic boundary conditions. In that case, sums of the 
form 5^(jj)( a i + a j) can be replaced by 4 £\ a i: and 

M" / , x M 

I ! 

,4 \A 



F 
~K 



' = \-f± A 1 1 |2 1 

1 V*cw + 7 A? (0)t c 2 0l ^ + £ 2(9.38) A 2 (0; 



E 



21^11^1 



Ai(0)t c oiAj(0)t c oj 



cos(0,: — 0A 



(29) 



Eq. (|2*9*|) is the most general form for the Ginzburg-Landau free energy functional considered 
in our calculations. In our simulations we allow both the amplitude \ip\ and the phase 9 of 
ip to undergo thermal fluctuations. 



B. Thermal averages 

As mentioned at the beginning of this section, at finite T we compute LDOS(u;,rj) by 
performing an average of LDOS(u;, {ipi}) over different configurations {ipi}. Those con- 
figurations are obtained assuming that the thermal fluctuations of {ipi} are governed by the 
Ginzburg-Landau (GL) free energy functional F described above. F is treated as an effective 
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classical Hamiltonian, and thermal averages (...) of quantities Q, such as LDOS(a;, ri), are 
obtained through 

{Q) = nfej^^iMi), (30, 

where Z is the canonical partition function, 

N 



\\d 2 ^e- F ' kBT . (31) 



i=i 



C. Estimate of the Kosterlitz-Thouless transition 

If amplitude fluctuations are neglected, the Hamiltonian ()29|) would correspond to an 
XY model on a square lattice. If the system is homogeneous, this XY model undergoes a 
Kosterlitz-Thouless transition at a temperature 

k B T c ~ 0.89J xy , (32) 
where Jxy is the coupling constant between spins: 

Hxy = -Jxy cos (^ " ( 33 ) 
From eqs. (|2*9^1 and ([3~3|). the XY coupling between sites i and j is given by 

■W*) = g^MI (34) 

\i(0)t c 0i\j(0)t c Qj 

If we approximate ipi(t) by the value that minimizes F' when fluctuations are neglected, 



\^{t)\ ~ >/9.38(l-*/*c0i)*c0i, (35) 

then 

, m 2(9.38) v /(l-t/t c0t )(l-t/t c0j ) 

J ^ (t) xmW) ' ( } 

which in the homogeneous case reduces to 

18.76(1 -t/t c0 ) 
•M*) ^ • (37) 



This result and eq. give 



Tc ~ ^— , (38) 

1 + Too/71 
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where 



Eq. (J38|) can also be rewritten as 



where 72 = Using A(0) = 1800 A and d = 10 A, we obtain j x = 172K. Finally, if 

choose Eq — 200 meV (for reasons given below), we obtain 72 = 13.54. 

Expressions (}3*8*j) and (}4T)j) will typically overestimate the phase-ordering (or Kosterlitz- 
Thouless) transition temperature T c . Both thermal fluctuations of and quenched disorder 
will generally reduce T c below these estimates. 



_ (0.89X18.76)^ 
71 " *»(0)* fl • (3J) 



tc * TTTV' (40) 

1 + t C 072 

' ~ "we 



D. Choice of parameters 

Next, we describe our choice of parameters entering both the microscopic model [Eq. (fT|)]. 
and that for thermal fluctuations [Eq. (|2*9*j)]. In a typical cuprate, such as underdoped 
Bi2212, the low-T superconducting gap is ~ 50 meV, the hopping integral t hop ~ 200 meV, 
A(0) ~ 1800A, and the pseudogap opens at T c0 ~ 200-K" ~ 20meV/kB- Also, the lattice 
constant of the CuC>2 lattice plane is a ~ 5.4A, while £ ~ 15A. If in eqs. (l26|) and (|27j) 
we choose E = th op = 200meV, then, using those expressions, we obtain |^(0)| = 0.25 
and t c o = 0.1. We can substitute these values into eq. (JHHJ) to obtain an estimate for the 
phase ordering temperature, namely T c = 130-fT. Our actual simulations, carried out in the 
presence of thermal fluctuations of the gap magnitude and quenched disorder, actually yield 
a lower T c , as expected. 

We have carried out calculations using this set of parameters, but also with smaller values 
of £0, in order to treat larger XY lattices. Suppose we wish to carry out a simulation on a 
16 x 16 XY lattice. If we use the parameters values described above, we would have a 48 x 48 
atomic lattice. To compute the density of states on this lattice, we would have to diagonalize 
4608 x 4608 matrices [see Eq. (|SJ)]. Each such diagonalization takes ~ 1 hour on a node for 
serial jobs of the OSC Pentium 4 Cluster, which has a 2.4 GHz Intel P4 Xeon processor. 
Since thermal averages require several hundred diagonalizations, a 16 x 16 XY lattice is too 
large using these parameters. If, however, we choose a smaller coherence length, we will 
have fewer atomic sites per XY cell, and hence a smaller matrix to diagonalize for a 16 x 16 
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XY lattice. In the BCS formalism, £o oc vf/\A\, where vf is the Fermi velocity. Thus, if £o 
is n times smaller than the experimental value, then, for fixed vf, A, and hence t C Q, will be 
n times larger than that value. 



E. Inhomogeneities 

As noted above, experiments show that in some cuprates the energy gap is spatially 



ovc, experiments snow 



inhomogeneous |2|, |3j, |4|, |5j, 13, |8j. Typically, in some spatial regions, which we call 
a-regions, the LDOS has a small gap and large coherence peaks, while in other regions, the 
/3-regions, the LDOS has a larger gap and reduced coherence peaks. The percentage of the 
area occupied by a and (3 regions, respectively, depends on the doping concentration. In 
BM2212, for example, a nearly optimally doped sample (hole dopant level ~ 0.18) has ~ 10 
% of the area occupied by (3 regions, while for an underdoped sample (hole dopant level 
~ 0.14), the areal fraction of the (3 regions is about ~ 50 % [4]. 

We introduce spatial inhomogeneities into our model by including a binary distribution of 
t c0 j's. Typically, we chose the smaller value of t c0i so that, for a homogeneous system, the gap 
Aj(0) resulting from our model [eq. EH approximately equals that observed in experiments 
(for further details, see discussion in the subsection entitled "choice of parameters"). We 
refer to XY cells with this small t c oi as a-cells. For the /3-cells, on the other hand, we 
assume a value t c0i K times large than that of the a-cells. We obtained our best results 
by choosing K = 3. We have carried out simulations considering both an ordered and a 
random distribution of /5-cells. 

To determine the distribution of Aj(0), we use the connection between the local superfluid 
density n Sji (T) and A f (T) implied by eqs. (fTSl). (ITT?|) and (1201: 

n s ,i(T) = mT)\ 2 = (9.38)16^™* (k^T^X^)' ^ 

Thus, at fixed but very low T, since |Aj(0)/fcsT c oi| 2 is independent of position according to 
our model [see Eq. (J22J], n s,i(T) oc 1/Af(0). Since the coherence peaks in the local density 
of states are observed to be lower where the gap is large, we will assume that t c0i and A^(0) 
are correlated according to the equation 

A, 2 (0) = ^t c0l (42) 
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where t c o and A 2 (0) are obtained from the observed bulk properties of the material un- 
der consideration. [For example, we typically obtain t c o from (|22|) where we take |Aj(0)| 
as the average of the low temperature gap observed in experiments, and A(0) = 1800A.] 
Substituting into eq. (J3B|) gives, for t « t c0i and t « t c0 j, 

JxY,i 3 oc — L= — (43) 
IV. COMPUTATIONAL METHOD 



A. Monte Carlo 

We compute thermal averages of several quantities, including LDOS(o;, T), using a 
Monte Carlo (MC) technique. Thus, we estimate integrals of the form (jHOj) using 

where N m is the number of configurations {ipi} used to compute the average, and the con- 
figurations {ipi} are obtained using the standard Metropolis algorithm |47|, |48j as we now 
describe. We first set the values of the t c oi and Aj(0) in each XY lattice cell as described 
in the previous section. This completely determines the GL free energy functional F [Eq. 
()29|).] We then set the initial values of ipi so as to minimize F. Next we perform attempts 
to change the value of each ipi by Si, where 5« is the complex number Si = S ijre + iS ijim , 
and 5i^ re and iS^im are random numbers with a uniform distribution in the range [— So, So}. 
We define a MC step as an attempt to change the value ipi on each of the XY cells. The 
value of Sq is in turn adjusted at each temperature so that attempts to change ip have a 
success rate of 50%. Attempts to change ipi are accepted with a probability exp(— AF/ k B T), 
where AF = F[ipi, ip 2 , . . . , ipi + Si, . . . , ip M ] -F[ipi,ip 2 , . . . ,ipi,. . . , iPm}- In this way, different 
configurations {ipi} are obtained. 

In order to select which of those confi gura tions {ipi} to use in (|44|). we first made an 
estimate the phase autocorrelation time r [28J , in units of MC steps, at each temperature. 
We chose r = min[r', 500], where r' is implicitly defined by 

<by = ? (45) 

14 



and c(r') is an space average of the phase autocorrelation function [48]: 

M 

c ( r ') = 7f E [(e^ {r,) e-^(°)) - (e^^Xe-^W)] . (46) 

Once we estimated r, we performed 20r MC steps to allow the system to equilibrate, then 
we carried out an additional lOOr MC steps at each T for each disorder realization. During 
those lOOr MC steps, we sampled {ipi} every r MC step, thus obtaining N m = 100 con- 
figurations to use in (jUJ) to estimate the quantities of interest. We also performed longer 
simulations, averaging over N m = 300 configurations to compute the LDOS, and N m = 5000 
configurations to compute 7, and the root-mean-square fluctuations [cr^] (defined be- 
low), obtaining virtually the same results as with N m = 100 configurations. 

When carrying out the simulation, we need a mapping between the sites of the XY lattice 
and those of the atomic lattice. To do this mapping, we divide the atomic lattice into regions 
of area £0 x £o- Each such region constitutes an XY cell. All atomic sites within such a cell 
are assigned the same value of the order parameter ipi. Clearly, the lattice constant £ of 
the XY lattice must be an integer multiple of the atomic lattice constant a . Thus, if our 
XY lattice has L 2 sites, then the atomic lattice has [L£ /ao] 2 sites. 

We diagonalize all matrices numerically using LAPACK 49] subroutine "zheev", which 
can find all of the eigenvectors and eigenvalues of a complex, hermitian matrix. We calculate 
the density of states by distributing the eigenvalues into bins of width Au. The delta function 
appearing in (fTB^) is approximated by 

6(x) = ~^— 2 , (47) 



where we choose e ~ Au ~ 0.0 It 



hop- 



B. Reducing finite size effects through inclusion of a magnetic field 

To reduce finite size effects on LDOS(co>, r), we use a method introduced by Assaad J^. 
The basic idea of this method is to break the translational invariance of th op through the 
substitution t hop — > Uj(L) in Eq. (j2J). This is done so as to improve convergence of the 
quantities of interest, such as LDOS(c<j,r), as a function of the size of the atomic lattice N 
|28j |. However, Uj(N) must still satisfy 

lim Uj(N) = -t^, (48) 

N— >oo 
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so that the original form of is recovered in the thermodynamic limit. 

Assaad showed that if one makes the substitution th. op — > Uj(N) through the inclusion 
of a finite magnetic field, the convergence of the density of states is greatly improved. The 
magnetic field enters through the Peierls phase factor: 



<^hop 6 11 1 



with 



2tt P 
$ 



Air) ■ dr. 



(49) 



(50) 



Here A[f) is the vector potential at r, $o — hc/e is the flux quantum corresponding to one 
electronic charge e, and the integral runs along the line from site i to site j. 

We use a gauge which allows periodic boundary conditions, and with which the flux 
through the atomic lattice can be chosen to be any integer multiple of $o 
i = (xe x a ,ye y a ), e x and e y are unit vectors in the x and y directions, and x and y are 
integers in the range [0, N — 1]. Then 



2irm 



0. 



if j = i± a e y , 
' "' !l- if J=i + a e x 
if j = i- a e x 
otherwise 



N 

27rm 
N y > 



and x — N — 1 , 
and x = 0, 



(51) 



where m is the number of flux quanta through the atomic lattice. We have chosen m = 1, 
so that the magnetic field in our system has the smallest non-zero value possible. 



V. RESULTS 



A. Zero temperature 

Fig. n shows the spatially averaged density of states, DOS(i^), obtained by summing 
the local density of states, LDOS(c<j,r), over all sites r on a 48 x 48 atomic lattice with 
homogeneous t c o at zero temperature. The zero temperature pairing strength is given by 
\i/j(0)\ = -\/9.38t c o, as shown by eqs. (jUJ) and (ffijj). For the case t c o = 0, the pairing strength 
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is zero, and we observe the standard Van Hove peak |52| for a two-dimensional tight-binding 
band at u — 0. For finite pairing strength we observe a suppression of the density of states 
near uj = 0, while strong coherence peaks occur at uj ~ |V>(0)|. 

In Fig. 2, we compare the density of states DOS(u;) for a 32 x 32 atomic lattice containing 
a single quantum of magnetic flux (q = 1), and a larger (48 x 48) atomic lattice containing 
no magnetic field (q = 0). 

Both systems are assumed homogeneous with t c o = 0.14. As can be seen, the two are 
very similar except at low \uj\, where the magnetic field is known to induce a change in the 
density of states Note also that the zero- field DOS(c<j) is less smooth than that of the 
lattice with one quantum of flux, even though the zero-field lattice is larger. In zero-field 
case we have determined the density of states using a bin width Au = 0.09, while in the 
finite-field case we used Au = 0.01. (The frequencies and widths are given in units of th op -) 
We have also carried out a similar calculation for q = 1 and a 48 x 48 atomic lattice; the 
results are similar to those shown for the 32 x 32 lattice except that the density of states at 
uj = is reduced by about a third. This Figure, and the results just mentioned, show that 
including the magnetic field is very useful in smoothing the density of states plots. 

Before presenting our results for inhomogeneous systems, we briefly describe our method 
of introducing inhomogeneities into our model. We work with atomic lattices of size L x L, 
in which the sites are divided into groups of 2 x 2. Each of these groups forms an XY cell, 
within which the superconducting order parameter ip is kept uniform. The value of i/j in 
each cell is determined by the GL free energy Eq. ()29j) . which in turn depends on the set of 
values {t c0i } and {A c oi}. Because t c0i and A c oi are correlated in our model, once we have the 
set {icoi}, the GL free energy is completely determined and ift at each cell can be computed 
through the MC method described above. 

In Fig. El (a) and (b), we show results for two inhomogeneous systems. Both systems 
consist of 48 x 48 atomic lattices in which a fraction cp = 0.11 of the XY cells are of the 
j3 type with t C Q = 0.42, while the remainder of the cells are of the a type, with t C Q = 0.14. 
The curves are spatial averages of the LDOS(o;,r) over the a and (3 cells. In (a), they 
correspond to a system in which the (3 cells form an ordered array , while the curves in part 
(b) correspond to a system in which the (3 cells are distributed randomly through the lattice. 
For comparison, Fig. |3fc) shows results of two homogeneous systems: one with all a cells 
and one with all (3 cells. 
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The dotted line in Fig. 01(a) represents an average of LD0S(a;,r) over the (3 cells. It 
differs significantly from the f3 curve of the homogeneous case, [dotted curve in Fig. Etc)]. 
Specifically, instead of the single sharp, and much higher, peak in the homogeneous (3 case, 
there is a lower peak which is shifted slightly to smaller \u\ and also has strong oscillations 
as a function of ui (probably because of the ordered arrangement of the (3 cells). The largest 
maximum of this oscillating peak is quite sharp, however, and occurs at a distinctly smaller 
energy than in the homogeneous case. 

The solid line in Fig. 0(a) corresponds to an average of the LDOS(u;,r) over a cells. It 
differs less from the homogeneous a system [solid curve in Fig. 0(c)] than in the (3 case: 
the main peak is not much shifted in energy, and it is slightly lower and broader than the 
homogeneous case. However, an additional peak does appear at the same position as the 
larger peak of the inhomogeneous (3 curve described above. 

In Fig. 0(b), we show the corresponding density of states plots for a system with randomly 
distributed (3 cells. In this case we observe that the LDOS(c<j,r), averaged over a cells, has 
slightly lower and broader peaks than that of the homogeneous a system shown in Fig. 0(c), 
but the peaks still occur at the same energy in both cases: u ~ 0.42. However, the average 
of the LDOS(u;, r) over the (3 cells is drastically different from the homogeneous (3 case: the 
main peak is greatly broadened, compared to the homogeneous (3 case. 

In Figs. E] and 01 we show representative ordered and disordered arrangements of a and 
f3 cells (for an 18 x 18 XY lattice), similar to those used in the calculations of Fig. 3(a) and 
3(b). In our density of states calculations for disordered arrangements, we typically average 
over about five realizations of the disorder, and use 24 x 24 XY lattices rather than the 
18 x 18 shown in the schematic picture. 

In Fig. 01 we show plots analogous to Fig. 01 but for a much larger concentration of f3 cells 
(cp = 0.89). Part (a) shows results for an ordered array of a cells immersed in a background 
of (3 cells. The simple, sharp peaks of the homogeneous (3 case [dotted curves in Fig. 0(c) 
and Fig. 0(c)] are split into two sharp peaks at a slightly smaller energy, while the sharp 
peaks of the homogeneous alpha regions, [solid curve in Fig. 0(c)] become even sharper and 
shifted toward higher energies, leading to a reduction in the density of states near u = 0. 
Also, in the inhomogeneous (3 curve of Fig. 0(a), a weak second peak appears at the same 
energy as one of the peaks in the inhomogeneous a curve. 

The case of a disordered distribution of a regions immersed in a background of (3 regions 
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is shown in Fig. |U] (b). The peaks of the curves corresponding to both the a and f3 regions 
become lower and broader than in the homogeneous cases, Fig. E^c). The peak in the (3 
curve occurs at approximately the same energy as in the homogeneous case. The corre- 
sponding inhomogeneous a peak, on the other hand, occurs at a higher energy relative to 
the homogeneous case. 

B. Finite temperatures 

We have carried out a fmite-T study for the system topology most similar to the exper- 
imental oiM'ljf : a random distribution of (3 regions immersed in a background of a regions. 
Calculated results for such a system at T = are shown in Fig. Efb). Because more matrix 
diagonalizations are required at finite T to obtain the relevant thermal averages, we work 
with 32 x 32 atomic lattices, instead of the 48 x 48 used at T = 0. Since the computational 
time needed for one diagonalization scales with the linear size L of the system like L 6 , each 
diagonalization takes about one-tenth the time in these smaller system. Fortunately, the 
reduction of finite-size effects achieved by introducing a magnetic field leads to good results 
even for this relatively small system size. This can be seen by comparing the t — results in 
Fig. [7J which are obtained for a 32 x 32 atomic lattice, to the corresponding results shown 
Fig. 0Kb) for a 48 x 48 atomic lattice. 

Besides the partial densities of states, we calculate several additional quantities at finite 
t: the effective superfluid density jit), the thermal- and space-averaged values of \ip\ in the 
a and f3 regions, and the relative fluctuations a\^\ of \ip\ averaged over each of those regions. 

We compute the superfluid density 7 by averaging the diagonal elements 7,™ (a = x, y) 
of the helicity modulus tensor 7. Thus, we compute 7 = {^ xx + 7 w )/2, where (37 1 

Ixx = Tj(^2( x i ~ %j) 2 JxY,ij c os(6>, - Oj)) - ~ Xj)J X Y,ij sin(0j - 9j)} 2 ) 

+Wt { ^ Xi ~ x ^ Jxy ^ sin( ^ - e ^ 2 - (52) 

Here Xi is the x coordinate of ith XY cell i, M is the total number of XY cells, JxY,ij is the 
effective XY coupling between XY cells and is given by Eq. 9i is the phase of ipi and () 
denotes a canonical average. 7 ra is defined by the analogous expression with X{ replaced by 
Hi. In our computations, we have set the lattice constant axY of the XY lattice to be unity. 
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The mean-square order parameter averaged over the a region is computed from 

KIVf>]« = ^EW 2 >> ( 53 ) 

where the sum is carried out over all M a XY cells of type a. [(IV'I 2 )]^ is defined similarly. The 
mean magnitude of the order parameter in the a and (3 regions, denoted [(iV'DJa an d [(|^|)]/3, 
are defined by an equation analogous to eq. ()53|). We compute the relative fluctuations [o"|^i] a 
of |^| within XY cells of type a from the definition 

(54) 

where the triangular brackets denote a thermodynamic average, and [...] denotes a space 
average over the a sites. is computed analogously. In systems with disorder, the 

square brackets denote a disorder average as well as a space average. 

Fig. [7| shows the partial LDOS(u;,r) averaged over a and (3 cells, at both t = and 
finite t. The systems shown have a fraction c@ — 0.1 of /3 sites randomly distributed. At 
t = the a regions show strong, sharp coherence peaks while the f3 regions have a larger 
gap but lower and broader peaks. When the temperature is increased to t — 0.015, the 
heights of both peaks are reduced, and their widths are increased, but the a peak is still 
quite sharp, because the system still has phase coherence. This temperature is still below 
the phase ordering temperature of t c ~ 0.03, as discussed below. As t is increased still 
further, to t — 0.035 and t = 0.055, the two density of states peaks broaden still further, 
there is scarcely any residue of a gap in the density of states, and there is now no sign of a 
real coherence peak in either the a or the (3 regions. 

In Fig.|Hl we show the superfluid density j(t), for the model just described but for various 
concentrations eg of the (randomly distributed) (3 cells. For eg = 0.1, the phase-ordering 
transition temperature t c ~ 0.03 in these units. Thus, of the plots in the previous Figure, 
two are below and two are above the phase-ordering transition. 

In Figs. El and we show the thermal, spatial, and disorder averages of over the a 
and (3 regions, denoted [(iV'DJa an d [(l^l)]/?, while Figs. [TTJ and H*2l show the corresponding 
averages of the root-mean-square fluctuations <j\m. A number of features deserve mention. 
First, the average is, of course, larger in the (3 regions than in the a regions, but the root- 
mean-square fluctuations are comparable in each of the two regions. Second, the increases in 
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the averages of |^| above the phase-ordering temperature is an artifact of a Ginzburg-Landau 
free energy functional as we now explain in detail. 

The asymptotic behavior of \ip\ 2 as t — > oo, for homogeneous systems, can be obtained in 
the following way : At very high temperatures t ^> t c0 , the first term in (J29j) goes like ~ \ip\ 2 , 
the second goes like ~ I^IV^ an d the third (coupling) term goes like ~ \ip\ 2 /t. We can then 
neglect the contribution of the third term, whence at large t the XY cells are effectively 
decoupled. The thermal average of \ijj\ 2 for an isolated cell is given by 

* = f °°w\m\ 2 eM-m 2 -gm , , 

] /;> W | e xp(-/H 2 -<#| 4 ) ' 1 ] 



In our case, 



/ = ip^ (56) 

and 

= Kl (iTi 

9 2(9.38)4A2(0)E t' 1 ' 

If / and g are real and positive, as in the present case, the integrals appearing in (J55|) can 
be carried out, with the result 

/i/|2\ / exp(-/ 2 /4#) 

(W> = -^ + ^m77v5)' () 

Here erfc(z) = 1 — erf(z), erf(z) being the gaussian error function. Using an asymptotic 
expansion 54j for erf(//2 v /^), applicable when / / \fg >> 1 as in the present case, we can 
show that 

lim(H 2 )-l//. (59) 

t— >oo 

Substituting / from eq. leads to 

]fa(\rP\ 2 )^ ^ X y )Eo -0.6 (60) 

where the last approximate equality is obtained using the parameters we have discussed 
above, namely K x ~ 2866 eV-A 2 , A(0) = 1800 A, t cQ = 0.14, E Q = 200meV. 
On the other hand, using eq. (J35j) . we obtain 



hm(|^| 2 )~0.2 (61) 

Thus, our model introduces an unphysical finite value of (\ip\ 2 ) at large t. This behavior has 
been observed in other studies of similar models while in other investigations this feature 
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is less obvious because of the parameters used|36j]. In our case, since we are interested in 
temperatures t < t c < t c0 or t ~ t c < t cQ , this unphysical high-temperature behavior should 
not be relevant to our calculations. 

Our results for [(IV'DL an d [(|^|)]^ suggest an explanation for one feature in the plots 
of the LDOS (see Fig. 7). Namely, the f3 peak generally occurs at higher uo than it would 
in a superconductor made entirely of (3 material. This shift occurs because, when the f3 
and a regions are mixed, [(IV'DJa is larger than its value in a homogeneous (3 system (as we 
further discuss below). This behavior of [(I^DU can be seen in Fig. [TDl where this quantity 
is plotted for different values of cp. Clearly, at low t, increases as c@ decreases. For the 
homogeneous (3 system, \ijj(t = 0)| = 1.29, as can be obtained directly from eq. this 
value is shown as an open circle at t — 0. This upward shift in the [(IV'I)]^ would be difficult 
to measure, since a pure (3 material may not exist. 

The behavior of [(I^DU has an analog, in our model, in the corresponding behavior in 
the a cells. Specifically, if a cells are the minority component in (3 host, tends to be 
substantially smaller than in pure a systems: the smaller the concentration c a = l — c@, the 
smaller the value of in those regions [see Fig. Ej. The behavior of in both a and 
(3 regions basically follows from our earlier discussion, according to which \ip\ 2 is larger in 
regions with a small gap. 

VI. DISCUSSION 

We have presented a phenomenological model for the temperature-dependent single- 
particle density of states in a BCS superconductor with a d x 2_ y 2 order parameter. Our 
model includes both inhomogeneities in the gap magnitude and fluctuations in the phase 
and amplitude of the gap. While some of these features have been included in previous 
models for the density of states (e. g. phase fluctuations in a homogeneous d-wave super- 
conductor, inhomogeneities in the gap magnitude at T = 0), our model is more general, and 
thus potentially more realistic for some cuprate superconductors. 

Our main goal is to examine the properties of an inhomogeneous superconductor, includ- 
ing many effects which are likely to be significant in real cuprate materials. The amplitude 
and phase fluctuations are treated by a discretized Ginzburg-Landau free energy functional, 
while the density of states is obtained by solving the Bogoliubov-de Gennes equations for a 
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superconductor with a tight-binding density of states and a d x 2_ y 2 energy gap. 

In all our calculations, we have assumed that the superconductor has two types of regions: 
a, with a small gap but a high superffuid density; and (3, with a large gap and small superffuid 
density. This assumption appears consistent with many experiments on the high-T c cuprates, 
especially in the underdoped regime If we assume that the minority component is of 
type P, embedded in an a host, we find that the local density of states at T = at the a 
sites has sharp coherence peaks, whereas that of the f3 sites is substantially broadened. This 
behavior is similar to experiment 

This description applies to a disordered distribution of (3 sites in an a host. If the (3 
sites are, instead, arranged on a lattice, the local density of states on the (3 sites is sharper, 
but also has distinct oscillations as a function of energy. Since such oscillations are absent 
in experiments, the actual (3 regions, if they exist as a minority component, are probably 
distributed randomly. 

In the reverse case of a regions embedded randomly in a f3 host, neither component has 
an extremely sharp density of states peak. While the a peak is still quite sharp, it is broader 
than the a peak in the /^-minority case. This result suggests that, if one component occurs 
only as isolated regions, its minority status tends to broaden its coherence peaks. 

Our results also show that the local density of states is strongly affected by phase fluctu- 
ations. This feature has already been found for a homogeneous d-wave superconductor j^, 
but here we demonstrate it in an inhomogeneous superconductor. The most striking effect of 
finite T is that the coherence peak in the a component disappears above the phase-ordering 
transition temperature T c . The (3 component does not show a coherence peak even at very 
low temperatures, but nonetheless this peak too is significantly broadened above T c . For T 
well above T c , there is no appreciable gap in the local density of states either at the a or 
the (3 sites. 

Our calculations include thermal fluctuations in the amplitude as well as the phase of 
tp. In general, thermal amplitude fluctuations seem to have only a minor influence on the 
local density of states. By contrast, the variations in |-0 1 due to quenched disorder (i. e., the 
presence of a and (3 regions in our model) strongly affect the local density of states, as we 
have already described. To check on the influence of purely thermal amplitude fluctuations, 
we have calculated the density of states of a homogeneous a superconductor with both phase 
and amplitude fluctuations, and have compared this to a similar calculation with only phase 



23 



fluctuations. We find that the additional presence of amplitude fluctuations has little effect 
on the density of states. 

To smooth the local density of states, we include in our density of states calculations 
a magnetic field equal to a flux hc/e in the entire sample area, following the method of 
Assaad [50]. This field greatly smooths the local density of states, which otherwise varies 
extremely sharply with energy, because of the many degenerate states of a finite sample at 
zero field. Our calculated density of states does, of course, correspond to a physical magnetic 
field, and thus differs slightly from that at zero field. For example, in a homogeneous system 
with a finite c?-wave gap, the LDOS(lg>) goes to zero as \u\ — > 0. By contrast, at finite field, 
the LDOS approaches a constant value at low |u;|. With no gap, our calculated DOS with 
nonzero field is indistinguishable from that of a conventional 2D tight-binding band (see 
Fig. 1), because the field is low (typically around 0.002 flux quanta per atomic unit cell). 
We conclude that the weak magnetic field very effectively smooths the calculated LDOS in 
a finite two-dimensional sample with a d-wave gap, but produces a density of states similar 
to that at zero field, except at very low \u\. 

Although we have included this magnetic field in the Bogoliubov-de Gennes equations, 
we have omitted it from the Ginzburg-Landau free energy functional, which is thus that 
of a zero-field system. As we now discuss, we believe that this numerical scheme should 
indeed converge to the correct physical result for zero magnetic field in the limit of a large 
computational sample. 

As noted earlier, we introduce the vector potential into the LDOS calculation in order 
to smooth the resulting density of states. In the limit of a large system, the effect of the 
vector potential, corresponding to a single quantum of flux, should become negligible, since 
the flux density becomes very small. This is already suggested by our calculated results for 
the two system sizes we consider (see Fig. 2 and the corresponding discussion). Even with a 
finite superconducting gap, the vector potential affects the LDOS very little, except at low 
energies; moreover, even this effect becomes smaller as the sample size increases. Therefore, 
in the limit of a large enough sample, our approach should give a very similar LDOS to 
one calculated with no vector potential. Hence, it is reasonable to use this approach in 
combination with a zero-field Ginzburg-Landau free energy functional to calculate the LDOS 
at finite temperatures. 

If we were to introduce a similar field into the Ginzburg-Landau functional, we believe 
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that it would have a substantial effect, not due to smoothing, on the phase ordering. There 
would be, not only the XY-like phase transition, as at zero field, but also additional phase 
fluctuations arising from the extra field-induced vortex. Since this extra vortex is absent at 
zero field, these effects would be irrelevant to the zero-field system we wish to model. By 
contrast, introducing a field into the Bogoliubov-de Gennes equations, as we do, provides 
desirable smoothing with little change in the calculated LDOS; moreover, even this slight 
change decreases with increasing sample size. Therefore, we believe that the best way to 
obtain a smooth LDOS at both zero and finite temperatures is to introduce the vector 
potential into the Bogoliubov-de Gennes equations, for smoothing purposes, but not to 
include it in the Ginzburg-Landau free energy functional. Our numerical results suggest 
that this procedure is indeed justified. 

One feature of our numerical results may seem counterintuitive. In our model, the a 
component is assumed to have a gap three times smaller than that of (3, but has a larger 
local superfluid density, i. e., a smaller penetration depth. We then find that the gap in the 
a region is smaller in a two-component system with both a and (3 regions, than it is in a 
pure a system. This counterintuitive result, however, emerges naturally from our discrete 
Ginzburg-Landau model, which is minimized if the quantities A i /(X i (0)T cOi ) are equal. For 
our model, Aj(0) is smaller in the small-gap material. It would be of interest if experimental 
evidence of this behavior were found in a real material. 

In our calculation, the LDOS is obtained from a Bogoliubov-de Gennes Hamiltonian 
whose parameters are determined by the Ginzburg-Landau functional. In fact, it should be 
possible to proceed in the opposite direction, and obtain the parameters of the functional 
from the LDOS. Specifically, the energy required to change a phase difference by a given 
amount depends on an integral over the LDOS. Thus, the calculation we have presented in 
this paper can, in principle, be made fully self-consistent. 
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FIG. 1: Zero temperature density of states DOS(cj) versus energy to for three homogeneous systems 
described by mean-field transition temperatures t C Q = 0, t C Q = 0.14, and t C Q = 0.42. Simulations 
were carried out on 48 x 48 atomic lattices, with a magnetic field included, as described in the 
text, to reduce finite size effects. Energies u are given in units of th op . 
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FIG. 2: Comparison of the T = DOS of a small system (32 x 32 atomic lattice) containing 
one quantum of magnetic field (q = 1) to that of a larger (48 x 48) system with no magnetic field 
(q = 0). Both systems are homogeneous with i c o = 0.14. Except at very low energies, the magnetic 
field produces little change in the shape of the curve. 
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FIG. 3: Comparison of the zero temperature DOS of mixed a- (3 systems (a) and (b), and pure 
systems (c). In the mixed systems, eg = 0.11 is the concentration of (3 cells (t c o = 0.42), immersed 
in a background of a cells (t c o = 0.14). (a) Ordered array of (3 cells; see Fig|lJ (b) Disordered 
configuration of (3 cells; see FigEl Curves are obtained by space- averaging LDOS(cj, r) over atomic 
sites within a or (3 cells, respectively. In the disordered case, averages were also carried out over 
five different realizations of the disorder, (c) DOS(w) for two pure systems containing only a and 
only (3 cells. We use a 48 x 48 atomic lattice; the size of an XY cell (a or (3) is 2 x 2. 
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FIG. 4: An 18 x 18 XY lattice in which the f3 cells form an ordered array: f3 cells (white squares) 
are immersed in a background of a cells (gray squares). Each XY cell corresponds to an area of 
x (o, and contains four atomic sites. 












































































































































































































































































































































































































































































































































































































































































FIG. 5: An 18 x 18 XY lattice with a disordered arrangement of (5 cells (white squares) immersed 
in a background of a cells (gray squares). Each cell corresponds to an area of £o x £o, and contains 
four atomic sites. This Figure contains a particular realization of disorder. Density of states results 
for disordered systems are averaged over five different disorder realizations. 
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FIG. 6: Same as Fig. 01 but with a high concentration cp = 0.89 of j3 cells in the mixed systems. 
The ordered configuration corresponds to an ordered arrangement of a cells within (3. 
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FIG. 7: Spatial average, at several different temperatures t, of the local density of states 
LDOS(u;,r) over two different types of cells: the a-cells, where t c o = 0.14, and the /3-cells, where 
t C Q = 0.42. The /3-cells occupy 10% (eg = 0.1) of the total area, while the a-cells occupy the rest. 
The simulations were performed using a 32 x 32 atomic lattice; the XY cells are 2x2 atomic 
cells. The phase ordering temperature for this system is t c ~ 0.03 (see the curve corresponding to 
op = 0.1 in Fig.©. 
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FIG. 8: Superfluid density *y(t) versus temperature t, for systems with different concentrations 
op of (5 cells distributed randomly over the atomic lattice. (3 cells have t C Q = 0.42, whereas a cells 
have t C Q = 0.14, However, the coupling constant between two nearest neighbor-cells (ij) includes, 
at low t, a factor i/v^cOi t C Qj, which results in a suppression of the superfluid density in systems 
with large concentrations of (5 cells [see Eq. (@HJ)]. 
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FIG. 9: Space, thermally, and disorder-averaged [(I^Dja; averaged over a cells, for systems with 
different concentrations cp of (3 cells, as described in the caption of Fig. |SJ In an a cell t C Q = 0.14 
while t C Q = 0.42 in a (3 cell. 
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FIG. 10: Same as Fig. |§1 but averaged over the (3 cells. 
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FIG. 11: Relative thermal fluctuations [<t\^\ (t)] a of \ip\, averaged over the a cells, for systems with 
different concentrations c@ of (5 as shown in Fig. |H1 
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FIG. 12: Same as Fig. ^2 but averaged over the j3 cells. 
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